////////////////////////////TCF/////////////////////////////
// TPC Coherent Fit (TCF)                                 //
// Author: Xianguo Lu, lu.xianguo@gmail.com               //
// Reference: Xianguo Lu, "Exploring the performance      //
// limits of the ALICE Time Projection Chamber and        //
// Transition Radiation Detector for measuring identified //
// hadron production at the LHC", dissertation at the     //
// University of Heidelberg (2013).                       //
////////////////////////////////////////////////////////////

#include "dEdxGausN.h"
#include "MN2D.h"

TTree *ftree = 0x0;
Int_t ix = -999;
Double_t rawp00=-999, rawp10=-999, rawp20=-999, rawp30=-999;
Double_t rawp01=-999, rawp11=-999, rawp21=-999, rawp31=-999;
Double_t rawp02=-999, rawp12=-999, rawp22=-999, rawp32=-999;

void FitBB(Double_t bbpar[], Double_t &realRes)
{
  TList *lgr=(TList*)gDirectory->Get("lgr");
  if(lgr){
    TH1D * hfit = (TH1D*) lgr->FindObject(Form("alephFITRev%d",dEdxGausN::LastfithIter()));
    if(!hfit){
      printf("mn2d FitBB no hfit!!\n"); lgr->ls(); exit(1);
    }
    
    TString omtit(hfit->GetTitle());
    for(Int_t iom=0; iom<dEdxGausN::NparBB(); iom++){
      const TString tmpf = omtit(0,omtit.First(',')-1);
      bbpar[iom]=tmpf.Atof();
      printf("mn2d FitBB fithRev mean used %s: %d %e\n", hfit->GetTitle(), iom, bbpar[iom]);
      omtit=omtit(omtit.First(',')+1,omtit.Length());
    }  
    const TString tmpf = omtit(0,omtit.First(',')-1);
    realRes = tmpf.Atof();
    printf("mn2d FitBB fithRev real Resolution used %s: %e\n", hfit->GetTitle(), realRes);
  }
  else{
    printf("mn2d FitBB no lgr! use default!\n"); //gDirectory->ls(); exit(1);
    memcpy(bbpar, dEdxGausN::IniBBpar(), sizeof(Double_t)* dEdxGausN::NparBB());
  }
}

void GetBinpar(Int_t &nbin, Int_t *binx, Double_t *binpar)
{
  printf("\nmn2d::GetBinpar ...\n");

  Int_t nx = 0;

  for(Int_t ii=0; ii<ftree->GetEntries(); ii++){
    ftree->GetEntry(ii);

    binx[nx] = ix;
    printf("mn2d::GetBinpar binx[nx] binx[%d], ix %d\n", nx, ix);

    //important to use rawp (natural ordering in fit), instead of p (particle ordering)
    //can be check by predraw/hnor to see whether everywhere = 1
    const Double_t rawps0[]={rawp00,rawp10,rawp20,rawp30};
    const Double_t rawps1[]={rawp01,rawp11,rawp21,rawp31};
    const Double_t rawps2[]={rawp02,rawp12,rawp22,rawp32};
    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      binpar[nx*MN2D::NparBin() + itype*dEdxGausN::NparType() + 0]=  rawps0[itype];
      binpar[nx*MN2D::NparBin() + itype*dEdxGausN::NparType() + 1]=  rawps1[itype];
      binpar[nx*MN2D::NparBin() + itype*dEdxGausN::NparType() + 2]=  rawps2[itype];
    }
    nx++;

    if(nx>nbin){
      printf("nx > nbin! %d %d\n", nx, nbin); exit(1);
    }
  }

  nbin = nx;
}


int main(int argc, char * argv[])
{
  TStopwatch localtime;
  localtime.Start(kTRUE);

  for(int ii = 0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }

  if(argc!=2){
    printf("argc!=4\n");return 1;
  }

  dEdxGausN::Ini();

  TString fin(argv[1]);
  //--------------------------

  TFile::Open(fin);

  //get preID from fith
  TH2D *hpreID = (TH2D*) gDirectory->Get("hpreID");

  //can NOT be done more than once! ==>
  ftree = (TTree*) gDirectory->Get("fith");
  if(!ftree){
    printf("mn2d no tree fith, use mn2d!!!!\n"); 
    ftree = (TTree*) gDirectory->Get("mn2d");
    if(!ftree){
      printf("mn2d no tree at all!\n");  gDirectory->ls(); exit(1);
    }
  }
  ftree->SetBranchAddress("ix",&ix);

  ftree->SetBranchAddress("rawp00", &rawp00);
  ftree->SetBranchAddress("rawp10", &rawp10);
  ftree->SetBranchAddress("rawp20", &rawp20);
  ftree->SetBranchAddress("rawp30", &rawp30);

  ftree->SetBranchAddress("rawp01", &rawp01);
  ftree->SetBranchAddress("rawp11", &rawp11);
  ftree->SetBranchAddress("rawp21", &rawp21);
  ftree->SetBranchAddress("rawp31", &rawp31);

  ftree->SetBranchAddress("rawp02", &rawp02);
  ftree->SetBranchAddress("rawp12", &rawp12);
  ftree->SetBranchAddress("rawp22", &rawp22);
  ftree->SetBranchAddress("rawp32", &rawp32);
 
  //use normalized 2dhist to deplete the physics influence between different bins
  const TString horiName("horiginal");
  printf("mn2d horiginal name: %s\n", horiName.Data());
  TH2D *horiginal=(TH2D*)gDirectory->Get(horiName);
  if(!horiginal){
    printf("no hh!\n");gDirectory->ls();exit(1);
  }

  //=======================================================
  //=======================================================
  Double_t bbpar[50];
  Double_t realRes = -999;
  FitBB(bbpar, realRes);

  //=======================================================
  //=======================================================

  Int_t tmpnx = 1000;
  Int_t *tmpbinx=new Int_t[tmpnx];
  Double_t *tmpbinpar=new Double_t[MN2D::NparBin()*tmpnx];
  GetBinpar(tmpnx, tmpbinx, tmpbinpar);

  const Int_t nx = tmpnx;
  Int_t binx[nx];
  Double_t binpar[MN2D::NparBin()*nx];
  memcpy(binx, tmpbinx, sizeof(Int_t)*nx);
  memcpy(binpar, tmpbinpar, sizeof(Double_t)*MN2D::NparBin()*nx);

  delete[] tmpbinx;
  delete[] tmpbinpar;

  printf("\nmn2d::GetBinpar done! %d\n", nx);
  //=======================================================
  //=======================================================
  MN2D::Ini(hpreID, horiginal, nx, binx, bbpar, binpar, realRes);

  //output: outmn2d_pre_draw2d.root
  MN2D::Draw2d("outmn2d_pre");

  //output: outmn2d_post_iBeta$(id)_save.root
  //id=100 is the final guided convergence results
  //id=0 is without regularization
  MN2D::Fit();

  //output: outmn2d_post_draw2d.root
  MN2D::Draw2d("outmn2d_post");
  
  MN2D::Clear();

  printf("mn2d done!!\n");  

  localtime.Stop();
  localtime.Print("m");
  return 0;
}
